#include <cassert>
#include <iostream>
#include <cstdlib>
#include <math.h>
#include "matriz.h"
using namespace std;

/***********************************
 * CONSTRUCTORES ******************
 *********************************/

Matriz::Matriz(){
    this->matriz = NULL;
    this->dim = 0;
	this->n = 0;
    this->m = 0;
}

Matriz::Matriz(unsigned int n, unsigned int m){
    assert((m > 0)&&(n > 0));
    this->dim = n*m;
    this->n = n;
    this->m = m;
    this->matriz = new double[this->dim];
    //inicializar en 0.
    for(unsigned int i=0;i<this->dim;i++)
        this->matriz[i] = 0;

}

Matriz::Matriz(unsigned int n, unsigned int m, unsigned int a){
    assert((m > 0)&&(n > 0));
    this->dim = n*m;
    this->n = n;
    this->m = m;
    this->matriz = new double[this->dim];
    //Adentro de cada lugar de la matriz va a ver un array de tamaño a, unicializado en 0
    double* datos = new double[a];
    for(unsigned int i=0;i<a;i++){
    datos[i] = 0;}
    
    for(unsigned int i=0;i<this->dim;i++)
        this->matriz[i] = *datos;

}



Matriz::Matriz(const Matriz& otra){
    if(this != &otra){
        this->dim = otra.dim;
		this->n = otra.n;
	    this->m = otra.m;
        this->matriz = new double[this->dim];
        for(unsigned int x = 0; x < this->dim; x++){
            this->matriz[x] = otra.matriz[x];
        }
    }
}

Matriz::~Matriz(){
    if(this->matriz != NULL)
        //delete [] this->matriz;
    this->matriz = NULL;
}

/***********************************
* OPERADORES         ***************
************************************/

/** deja en este parametro la suma de este con el otro */
Matriz Matriz::suma(const Matriz& otra){
    assert((this->n==otra.n)&&(this->m==otra.m));
	Matriz res(this->n, this->m);
    for(unsigned int i=0;i<this->n;i++)
        for(unsigned int j=0;j<this->m;j++)
                res.setPos(i,j, (*this)(i,j)+otra(i,j) );

	return res;
}

/** deja en este parametro la resta de este con el otro */
Matriz Matriz::resta(const Matriz& otra){
    assert((this->n==otra.n)&&(this->m==otra.m));
	Matriz res(this->n, this->m);
    for(unsigned int i=0;i<this->n;i++)
        for(unsigned int j=0;j<this->m;j++)
                res.setPos(i,j, (*this)(i,j)-otra(i,j));

	return res;
}

/** devuelve el vector de largo dim del resultado de la multiplicacion*/
Matriz Matriz::producto(const Matriz& otra){
    assert(this->m==otra.n);
	Matriz res(this->n, otra.m);
    for(unsigned int i=0;i<(this->n);i++) {
        //guardo la fila nueva tmpmente para no sobreescribri nada
        //double aux[this->m];
        for(unsigned int j=0;j<(otra.m);j++) {
                double val = 0;
                for(unsigned int k=0;k<(this->m);k++)
                        val += ( (*this)(i,k) * otra(k,j) );
                //aux[j] = val;
				res.setPos(i,j,val);
        }
        /*for(unsigned int j=0;j<this->m;j++)
                this->setPos(i,j,aux[j]);
		*/
    }
	
	return res;
}

/** deja en esta matriz el resultado del producto por un escalar*/
Matriz Matriz::producto(const double lambda){
	Matriz res(this->n, this->m);
    for(unsigned int i=0;i<this->n;i++)
        for(unsigned int j=0;j<this->m;j++)
                res.setPos(i,j, (*this)(i,j)*lambda);

	return res;
}


Matriz Matriz::triangular(){
	Matriz res = (*this);

    unsigned int i=0,j=0;
    //1)Triangulamos la matriz junto con b
    while(i<res.n && j<res.m) {
        //1.1)Hacemos el pivoteo parcial
        unsigned int maxRow = i;
        for(unsigned int k=i+1;k<res.n;k++) {
            if(fabs(res(k,j))>fabs(res(maxRow,j))) {
                maxRow=k;
            }
        }
        //if(res(maxRow,j)==0) return false;
        //1.2)Cambiamos la fila por la de maximo valor en la columna actual
        res.swapRows(i,maxRow);
        //1.3)Hacemos el siguiente paso de la triangulacion
        for(unsigned int k=i+1;k<res.n;k++) {
            double valAux=res(k,j)/res(i,j);
            for(unsigned int l=j;l<res.m;l++) {
                if(l==j) {
                    res.setPos(k,l,0);
                } else
                    res.setPos(k,l,res(k,l)-valAux*res(i,l));
            }
        }
        i++;
        j++;
    }
	return res;
}

bool Matriz::getLU(Matriz& L, Matriz& U){
	float u,l;
	Matriz A = (*this);
	if(A(1,1) == 0)
		return false;

	L.setPos(1,1,1);
	U.setPos(1,1,A(1,1));

	for(unsigned int r = 2; r<=this->n; r++){
		U.setPos(1,r,(A(1,r)/L(1,1)));
		L.setPos(r,1,(A(r,1)/U(1,1)));
	}

	for(unsigned int i = 2; i<=(this->n -1); i++){
		L.setPos(i,i,1);
		u = A(i,i);
		for(unsigned int k = 1; k<i-1;k++){
			u -= L(i,k)*U(k,i);
		}
		if(u == 0)
			return false;

		U.setPos(i,i,u);

		for(unsigned int j = i + 1; j<=this->n; j++){
			u = A(i,j);
			for(unsigned int k = 1; k<=i-1;k++){
				u -= L(i,k)*U(k,j);
			}
			u = (1/L(i,j))*u;
			U.setPos(i,j,u);
			
			l = A(i,j);
			for(unsigned int k = 1; k<=i-1;k++){
				l -= L(j,k)*U(k,i);
			}
			l = (1/L(j,i))*l;
			L.setPos(j,i,l);
		}
	}
	
	L.setPos(this->n,this->n,1);
	u = A(this->n,this->n);
	for(unsigned int r = 1; r<=this->n-1; r++){
		u -= L(this->n,r)*U(r,this->n);
	}
	U.setPos(1,1,A(1,1));

	return true;
}

Matriz Matriz::triangular(double* be){
	//se copian los valores para modificarse
    Matriz res = (*this);

    double b[this->n];
    for(unsigned int i=0;i<this->n;i++)
        b[i]=be[i];
    unsigned int i=0,j=0;
    //1)Triangulamos la matriz junto con b
    while(i<res.n && j<res.m) {
        //1.1)Hacemos el pivoteo parcial
        unsigned int maxRow = i;
        for(unsigned int k=i+1;k<res.n;k++) {
            if(fabs(res(k,j))>fabs(res(maxRow,j))) {
                maxRow=k;
            }
        }
        //if(res(maxRow,j)==0) return false;
        //1.2)Cambiamos la fila por la de maximo valor en la columna actual
        res.swapRows(i,maxRow);
        double tmp=b[i];
        b[i]=b[maxRow];
        b[maxRow]=tmp;
        //1.3)Hacemos el siguiente paso de la triangulacion
        for(unsigned int k=i+1;k<res.n;k++) {
            double valAux=res(k,j)/res(i,j);
            for(unsigned int l=j;l<res.m;l++) {
                if(l==j) {
                    res.setPos(k,l,0);
                } else
                    res.setPos(k,l,res(k,l)-valAux*res(i,l));
            }
            b[k]=b[k]-valAux*b[i]; //tmb en b!
        }
        i++;
        j++;
    }
	for(unsigned int i=0;i<this->n;i++)
        be[i]=b[i];
	//(*be) = b;
	return res;
}

bool Matriz::resolver(double* be, double* x) {
    //se copian los valores para modificarse
    Matriz aux = this->triangular(be);

    //Resolvemos la matriz ya triangulada
    for(int k=aux.n-1;k>=0;--k) {
        double val=be[k];
        for(unsigned int l=k+1;l<aux.m;l++) {
            val-= ( aux(k,l)*x[l] );
        }
        x[k]=val/aux(k,k);
    }
    //esto de aca de va a borrar
    /*cout << "LA MATRIZ QUEDA COMO: " << endl;
    for(unsigned int k=0;k<aux.size();k++) {
        for(unsigned int l=0;l<aux.size();l++) {
            cout << aux(k,l) << '\t';
        }
        cout << endl;
    }
    cout << "Y LA SOLUCION ES: [";
    for(unsigned int k=0;k<aux.size();k++)
        cout << x[k] << ", ";
    cout << "]" << endl;*/
    return true;
}

//Matriz operator+(const Matriz& a, const Matriz& b){
//    assert(a.n == b.n);
//    assert(a.m == b.m);
//    Matriz res = a.suma(b);
//    
//    return res;
//}
//
//Matriz operator-(const Matriz& a, const Matriz& b){
//    assert(a.n == b.n);
//    assert(a.m == b.m);
//    Matriz res= a.resta(b);
//    
//    return res;
//}
//
//Matriz operator*(const double& lambda, const Matriz& a){
//    return a.producto(lambda);
//}
//
//
//Matriz operator*(const Matriz& a, const double& lambda){
//    return lambda*a;
//}
//
//Matriz operator*(const Matriz& a, const Matriz& b){
//	assert(a.m==b.n);
//    Matriz res= a.producto(b);
//    return res;
//}

ostream& operator<<(ostream& out, const Matriz& m) {
    for(unsigned int i=1;i<=m.n;i++) {
        for(unsigned int j=1;j<=m.m;j++) {
            out << m(i,j) << '\t';
        }
        out << endl;
    }
    return out;
}

/**********************************
* GETTERS    **********************
**********************************/
double Matriz::getId(unsigned int fila, unsigned int columna) const {
    return (fila-1)*this->m+(columna-1);
}

/** devuelve el valor en la posicion (fila,columna) */
double Matriz::operator()(unsigned int fila,unsigned int columna) const{
//cout<<fila<<"  -  "<<columna<<"  /  "<<this->n<<"  -  "<<this->m<<endl;
    assert(fila <= this->n && columna <= this->m);
    return this->matriz[(fila-1)*this->m + (columna-1)];
}

double Matriz::operator()(unsigned int id) const{
    return this->matriz[id];
}


/** devuelve el tamanio de la matriz */
unsigned int Matriz::size() const{
    return this->dim;
}

/**********************************
* SETTERS    **********************
**********************************/
/** setea el valor en la posicion (fila,columna) */
void Matriz::setPos(unsigned int fila, unsigned int columna, double valor) {
    assert(fila <= this->n && columna <= this->m);
    this->matriz[(fila-1)*this->m + (columna-1)] = valor;
}

void Matriz::setPosCampoN(unsigned int fila, unsigned int columna, double* valor) {
	assert(fila <= this->n && columna <= this->m);
    this->matriz[(fila-1)*this->m + (columna-1)] = *valor;
}


void Matriz::setPosInId(unsigned int id, double valor) {
        this->matriz[id]=valor;
}

/*********************************
* COMPARACION Y ASIGNACIONES *****
*********************************/
bool Matriz::operator==(const Matriz& other) const {
    if( this != &other) {
        if((this->n != other.n)||(this->m != other.m)) {
            return false;
        } else {
            bool igual = true;
            for(unsigned int i = 1; i <= this->n && igual; i++) {
                for(unsigned int j = 1; j <= this->m && igual; j++) {
                    igual = igual && (*this)(i,j) == other(i,j);
                }
            }
            return igual;
        }
    }
    return true;
}

bool Matriz::operator!=(const Matriz& other) const{
    return !((*this) == other);
}

Matriz& Matriz::operator=(const Matriz& other){
    if( this != &other) {
        this->dim = other.size();
		this->n = other.n;
		this->m = other.m;
        if( this->matriz != NULL)
            delete [] this->matriz;
        this->matriz = new double[this->n*this->m];
        for(unsigned int x = 1; x <= this->n*this->dim; x++)
            this->matriz[x] = other.matriz[x];
    }
    return *this;
}

void Matriz::swapRows(unsigned int i, unsigned int j) {
    assert(i<this->n && j<this->n);
    if(i!=j) {
        for(unsigned int k=1;k<=this->m;k++) {
            double tmp = (*this)(i,k);
            setPos(i,k,(*this)(j,k));
            setPos(j,k,tmp);
        }
    }
}

bool Matriz::SOR(int n, Matriz A, float* b, float* X, float w, float tol, int iterMax, float* res){
	int k = 1;
	//float* res = new float[n];
	float aux;
	while(k<iterMax){
		for(int i=1; i<n;i++){
			aux = w;
			for(int j=1;j<=i-1;j++){
				aux -= A(i-1,j-1)*res[j-1];
			}

			for(int j=i+1;j<=n;j++){
				aux -= A(i-1,j-1)*X[j-1] + b[i-1];
			}
			//aux += b[i];
			res[i-1] = (1-w)*X[i-1] + (aux/A(i-1,i-1));
		}
		
		aux = 0;
		for(int r = 0; r<n; r++){
			aux += (res[r] - X[r])*(res[r] - X[r]);
		}
		if(sqrt(aux)<tol){
			return true;
		}

		k++;
		for(int r = 0; r<n; r++){
			X[r] = res[r];
		}
	}

	cout<<"se paso de rosca"<<endl;
	return false;
}


